Interaction effects in transport through molecular monolayers 
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We theoretically investigate the effect of inter-molecular Coulomb interactions on transport 
through molecular monolayers (or other devices based on a large number of nanoscale conductors 
connected in parallel). Due to the interactions, the current through different molecules become cor- 
related, resulting in distinct features in the non-linear current-voltage characteristics, as we show by 
deriving and solving a type of modified master equation, suitable for describing transport through an 
infinite number of interacting conductors. Furthermore, if some of the molecules fail to bond to both 
electrodes, charge traps can be induced at high voltages and block transport through neighboring 
molecules, resulting in negative differential resistance. 

PACS numbers: 85.65.+h, 87.15.hg, 85.35.Gv, 85.35.-p, 



Introduction. It is well known that when N macro- 
scopic resistors, each with a conductance G, are con- 
nected in parallel, the total conductance is given by 
Gt = NG. However, this is no longer true if the currents 
flowing through the individual resistors are correlated, 
i.e., if the current through one resistor is affected by the 
current through its neighbors, e.g., due to Coulomb inter- 
action between charge carriers. Such correlation effects 
can be expected to be important in nanoscale systems, 
where the conductors are very close to each other and 
where, because of the quantized nature of the electron 
charge, current flow can be associated with significant 
fluctuations of the charge on the resistors. 

One interesting example is molecular electronic devices 
based on a self-assembled molecular monolayer, sand- 
wiched between metallic electrodes [IHH . Compared with 
single-molecule junctions, monolayer devices offer better 
reproducibility and stability, and the larger currents are 
easier to measure. They have been used to investigate a 
number of interesting molecular transport effects, such as 
negative differential resistance (NDR) jij], switching Q, 
and spin-selective tunneling 0, Q] . Although the exper- 
iments are usually interpreted within a single-molecule 
picture, direct experimental comparisons between mono- 
layer and single-molecule junctions Q have found large 
differences, i.e., Gt ^ NG. Such differences have been 
discussed in terms of static changes to the local molec- 
ular environment in a monolayer device, e.g., due to re- 
arrangement of molecular or surface charges, or interac 
tions between constant molecular dipolc moments 



In this work, we investigate the dynamic transport ef- 
fect resulting from Coulomb interactions between charges 
being transported through neighboring molecules in a 
monolayer. The inter-molecular Coulomb interactions 
not only lower the conductance (Gt < NG), but qual- 
itatively change the non-linear current-voltage charac- 
teristics of the device, as we show by deriving a type 
of modified master equation for the nonequilibrium cur- 
rent, as well as the voltage-dependent charging of the 
monolayer. If the source and drain tunnel couplings dif- 



fer, the inter-molecular Coulomb interactions give rise 
to a voltage-asymmetric current-voltage characteristic, 
I(V) 7^ — J(— V), the shape of which provides an ex- 
perimental fingerprint of the interactions. Furthermore, 
if some molecules form bonds only with one electrode, 
we show that NDR can occur, since charge traps are 
formed within the layer and, through the inter-molecular 
Coulomb interaction, block transport through neighbor- 
ing molecules. Knowledge of the generic transport signa- 
tures of inter-molecular Coulomb interactions, derived in 
the present work, should be very helpful when interpret- 
ing data from molecular monolayer devices and trying 
to separate the genuinely single-molecule transport ef- 
fects. For simplicity, we use terminology like "molecule" 
and " monolayer" , but the results are of much more gen- 
eral importance, and apply, for example, to transport 
through many quantum dots connected in parallel, car- 
bon nanotubes bundled together in a rope, or arrays of 
nanoparticles, just to give a few examples. The new type 
of master equation we develop, which describes the cor- 
related non-equilibrium charge distribution among inter- 
acting conductors, should form a useful starting point 
also for studies of a wider range of interacting mesoscopic 
systems. 

There has been a significant amount of prior work 
investigating the transport effects of an inter-molecular 



tunnel coupling, see e.g., Refs. I12H15I . However, in con 



trast to Coulomb interactions, tunnel couplings decay ex- 
ponentially with the molecular separation. Several works 
have also studied transport through two parallel-coupled 
quantum dots 

EMI, 

including the effects of inter-dot 

Coulomb interactions. 

Basic physical picture and model. Figure [T] shows a 
sketch of a molecular monolayer between metallic source 
and drain electrodes. In Fig. [Tf b) , we focus on two 
molecules within the layer and illustrate the basic trans- 
port effect of inter-molecular Coulomb interactions, and 
indicate the molecular transmission resonances (the ver- 
tical direction therefore indicates a direction both in real 
space and in energy space). When an electron tunnels 
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FIG. 1. (Color online), (a) Sketch of the molecular mono- 
layer device, (b) Zoom in on two of the molecules in the up- 
per panel, together with an energy-space representation of the 
molecular transmission resonances. Due to the intermolecu- 
lar Coulomb interaction, the resonance seen by an electron 
tunneling into a molecule is shifted by U if one of its neigh- 
bors is charged, (c) Level diagram showing tunneling into 
the LUMO of a single molecule. The LUMO is shifted by 
kll when k neighboring molecules are occupied (higher-lying 
dashed lines). 



into a molecule, it experiences the electric field from 
charges localized on neighboring molecules, which leads 
to a shift of the effective transport resonance. Therefore, 
when the bias voltage, V, is just barely large enough 
to allow resonant transport through a molecular orbital 
(the HOMO or LUMO is inside the bias window), trans- 
port can still be suppressed if neighboring molecules are 
charged. In the limit of weak tunnel coupling, where 
transport is dominated by single-electron tunneling, and 
at low temperatures, we therefore expect a single molec- 
ular orbital to give rise to multiple transport resonances: 
One when the bias is large enough to allow tunneling into 
the LUMO or out of the HOMO, and additional ones 
when this becomes possible also when k = 1, 2, . . . neigh- 
boring molecules are charged, see also Fig. QTc) . This is 
in constrast to the result without inter-molecular inter- 
actions, where only a single resonance occurs. Below we 
calculate the detailed current- voltage curves for different 
interaction strengths. 

The distance between molecules in a monolayer de- 
pends on the molecular species and contact materials, 
but is typically d ~ 1 nm in a dense monolayer. A 
rough estimate gives the interaction between electrons on 
neighboring molecules as U ~ e 2 / (4ndeoe r ) ~ 1 eV/e r . 
The relative permittivity, e r , within the monolayer may 
be rather large and the interaction strength further de- 
creased by screening from conduction electrons in the 
electrodes However, a U of the same order of mag- 
nitude as the thermal energy at room temperature still 
seems reasonable. In other mesoscopic systems, U is 
likely much smaller (U = 0.15 — 0.4 meV between strands 



of nanotubes within a nanotube rope was found in recent 
experiments [20I 21 1), but can still dominate the trans- 



port physics at low temperatures. 

To focus on the interaction effects, we study the sim- 
plest possible molecular model, including only a sin- 
gle spin-less orbital (chosen as the LUMO for definitc- 
ness). The Hamiltonian of the monolayer device is 
H = Ei(#ML + Ht) + H Iea , where 
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Here, is the energy of the LUMO of molecule i, which 
has occupation rii = d\di, and Uij is the Coulomb charg- 
ing energy between electrons on dots i and j (in all 
calculations we include only nearest neighbor interac- 
tions). The source (r = S) and drain (r = D) elec- 
trodes are described by H ves , where 

Upr — Cp r Cp r is the 
number operator for electrons in state p. The electrons 
in the electrodes are as usual for good metals modelled 
as being non-interacting. describes tunneling be- 

tween electrode r = S,D and molecule i, which takes 
place with amplitude t r i. We assume both the tunnel 
amplitude and the electrode densities of states, p ri to 
be energy-independent. In this case, the tunnel rates, 
r r j = 2irp r \t r i\ 2 /H, which set the inverse time-scale for 
single-electron tunneling, are also energy- independent. 

Transport through a homogeneous layer. We consider 
first a homogeneous monolayer, i.e., Uij — U, T r i = T r , 
and ti = e. We also focus on the regime where transport 
is dominated by single-electron tunneling, which is the 
case in the weak tunnel coupling regime, HT < fc^T, 
where T is the temperature. We arbitrarily pick one 
molecule in the monolayer and denote it by Ml. Our 
aim is now to calculate the stationary current flowing 
through the LUMO of Ml. In the single-electron tunnel- 
ing regime, the current is carried by processes in which 
a single electron tunnels either from one electrode into 
Ml, or from Ml into one electrode, with rates propor- 
tional to the probability that Ml is empty, P^, or full, 
PW, respectively. However, because of U, the rates also 
depend on the occupation of the TV nearest neighbors of 
Ml, which we denote M2. We therefore consider prob- 
ability distributions of the form ■ The probability 
that Ml is empty (m = 0) or occupied (m = 1), while 
k of its neighbors are occupied. The net current from 
reservoir r into Ml is then 
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lh)/h B T), /+0) is the 
f (a;), and p r = ±eV/2 
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FIG. 2. (Color online). I(V), (a), and dI(V)/dV, (b), 
for e = lOfcsT and increasing strength of the inter-molecular 
Coulomb interaction, U = 0, . . . ,AksT. We have assumed a 
square lattice configuration, where each molecule has N — 4 
nearest neighbors. / and dl/dV are given in arbitrary units 
(au) since I oc F in the single-electron tunneling regime. The 
inset in (b) shows the linear conductance, Gt, as a function 
of the position of the LUMO relative to the Fermi energy. 

for r = S/D is the chemical potential of electrode r. 
We assume that the bias drop occurs symmetrically at 
the two tunnel barriers, such that e is independent of V . 
The total current is / = NmId = —NmIs, where Nm is 
the total number of molecules in the monolayer. 

The remaining problem is to find a master equation to 
determine Pj. under non-equilibrium conditions, i.e., 
to find the voltage-induced charging of the interacting 
monolayer. In the supplementary online information [22| . 
we present the complete master equation and the full 
derivation, and also compare solutions within different 
approximation schemes. Here we merely mention that a 
master equation for Pj. involves tunneling into and out 
of both Ml (changing m) and M2 (changing k). The lat- 
ter terms are proportional to the conditional probabilities, 
Pj.™! ■ The probability that M2 is occupied (n = 1) or 
empty (n = 0) and has k' occupied neighbors, given that 
Ml is occupied (m = 1) or empty (m = 0) with k occu- 
pied neighbors. The master equation for PiTif 1 ^ in turn 

involves higher order conditional probabilities, P^"jj) ■ 
In the results presented here, we close the hierarchy of 
equations by neglecting explicit three-charge correlations 
(replacing Pj.™? 1 ^) —> P^'S/)- One is then left with solv- 
ing a non-linear master equation for Puur n \ which have 
to be inserted into the equation for P^ ■ 

Figure [2] shows the results of the master equation cal- 
culations for U = 0, . . . , 4k B T. When U = 0, I(V) shows 
a thermally broadened step when the LUMO enters into 
the bias window, which acquires some additional broad- 
ening for U < ksT. When U > ksT, several current 



FIG. 3. (Color online). |/|, (a), and dl/dV, (b), as a function 
of V (thick lines) or —V (thin lines), for e = 10fc B T, U = 
4fcsT, and N = 4. We vary the ratio between the source 
and drain tunnel couplings, a = Fs/Td, while keeping F = 
TsTd/{Ts + F_d) fixed, which fixes the value of the current 
plateau at large V. 



steps (multiple sidepeaks in dl/dV) can be discerned. 
They are separated in voltage by U and correspond to 
the condition that charging Ml becomes possible, also 
when k = 1, 2, . . . , N of M2 are charged. There are only 
N satellite conductance peaks, since the maximum inter- 
action cost of adding one electron is NU, and the peaks 
are always equidistant. This makes it possible to ex- 
perimentally distinguish sidepeaks originating from inter- 
molecular Coulomb interaction from similar features re- 
lated to higher lying orbitals or vibronic excitations. For 
all values of U , the current saturates at the same value, 
corresponding to a fully conducting monolayer. The in- 
set of Fig. [H[b) shows the linear conductance as a func- 
tion of the position of the LUMO relative to the elec- 
trode Fermi energy. Unless the LUMO is close to reso- 
nance, the charging of the monolayer is rather small and 
the inter-molecular Coulomb interactions are relatively 
unimportant. 

Figure |3] shows I(V) and dI(V)/dV for a monolayer 
with a larger coupling to the source than the drain, 
Ts > I\d. For {7^0, this introduces an asymmetry 
into the current voltage characteristics, I(V) ^ — /(— V). 
The reason is that there is now a larger charging of the 
monolayer for V > compared with V < since, for pos- 
itive bias, electrons can easily tunnel into the monolayer 
(from the source), but cannot easily escape again (to the 
drain), while the situation is reversed for negative bias. 
A larger charging increases the effects of interactions and 
causes the I(V) curve to become flatter and dl/dV to 
show more pronounced sidepeaks. For [7 = 0, asymmet- 
ric tunnel couplings do not introduce any asymmetry into 
the I(V) curve. There are several other possible reasons 
for an asymmetric /(V), but the rather special form seen 
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FIG. 4. (Color online). I(V) for e = 517/2, N = 4, and 
successively lowered T. One molecule (out of nine) has Fd — 
Fs/100. In the magenta curve, we have furthermore assumed 
the LUMO of the asymmetrically coupled molecule to be 8e = 
U/4 higher in energy compared with the others. The inset 
zooms in on the first current step. 



in Fig. [3] nonetheless provides a fingerprint of the inter- 
molecular Coulomb interaction. 

NDR in inhomogeneous layers. Finally, we investigate 
the effects of disorder within the molecular monolayer. 
Here, it is no longer possible to use our master equation 
approach, which assumes a homogeneous monolayer. In- 
stead, we take a finite number of molecules, which may 
all have different properties, and diagonalize the corre- 
sponding many-body Hamiltonian exactly. The many- 
body eigenstates can then be used in a standard mas- 
ter equation approach [23| to calculate the current. Due 
to the rapidly growing size of the Hilbert space, only a 
rather small number of molecules can be included. How- 
as we show in the supplementary online informa- 



22j for the case of a homogeneous monolayer, using 



ever, 
tion 

3x3 molecules and periodic boundary conditions gives 
excellent agreement with the infinite monolayer master 
equation. 

Figure [4] shows the result of these calculations, where 
we have assumed one molecule (out of nine) to be only 
weakly coupled to the drain electrode, Tjj = Ts/lOO. 
The asymmetrically tunnel coupled molecule now acts 
as a charge trap, since for V > 2e it can be filled 
with an electron tunneling in from the source, which is 
then prevented from escaping to the drain. The inter- 
molecular Coulomb interaction means that the occupa- 
tion of the charge trap prevents transport through neigh- 
boring molecules, unless the bias voltage is increased fur- 
ther. As the temperature is reduced, this leads to a 
weak NDR at the first current plateau. A much stronger 
NDR effect, which is repeated on every current plateau, 
is obtained if the LUMO of the asymmetrically coupled 
molecule lies at a slightly higher energy compared with 
the other molecules (see magenta curve in Fig. 2]). This 
is a realistic scenario since the asymmetrically coupled 
molecule experiences less screening from conduction elec- 
trons in the electrodes and, in general, a different chem- 
ical and electrostatic environment. The reason for the 
larger NDR is that, due to the misalignment, a large cur- 
rent is first allowed to flow through the symmetrically 



coupled molecules, before the charge trap becomes filled 
at higher voltages. For the configuration in Fig. [4] with a 
molecule weakly coupled to the drain, the I{V) curve will 
be asymmetric, with no NDR for negative bias. However, 
in a large monolayer, if approximately the same number 
of molecules are weakly coupled to the source and to the 
drain, the voltage symmetry is approximately restored 
with NDR for both bias polarities. 

Conclusions. We have studied electric transport 
through a large number of mesoscopic conductors con- 
tacted in parallel, and shown that Coulomb interactions 
between charge carriers on neighboring conductors give 
rise to distinct features in the non-linear current-voltage 
characteristics. Our detailed calculations allow these 
interaction-induced features to be identified and isolated 
from the single-device properties, which is essential when 
using electric transport as a spectroscopic tool. Further- 
more, we have shown that interactions can give rise to 
negative differential resistance if some of the conductors 
are significantly coupled to only one electrode. Our re- 
sults may aid in the interpretation of past and future ex- 
periments, for example on molecular monolayer devices. 
There has been much recent interest in molecular ther- 
moelectric devices 24 , 2{| because of the sharp nature of 
the transport resonances [26| . The results of our study, 
showing an effective interaction-induced broadening of 
the transport resonances, is likely to have a large impact 
on the thermoelectric efficiency. Finally, the general the- 
oretical framework, derived in the supplementary infor- 
should form a useful basis for future studies 
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of transport through interacting mesoscopic conductors. 
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SUPPLEMENTARY INFORMATION 
General master equations 

We start by deriving general expressions for the cur- 
rent through the interacting molecular monolayer, and its 
noncquilibrium charging, focusing on the single electron 
tunneling regime. We consider a homogeneous monolayer 
and (arbitrarily) pick one molecule (denoted Ml). We 
want to calculate the steady-state current flowing from 
electrode r into Ml, which is given by (see Eq. (4) in the 
main paper) 



N 

<E 

k=0 



W+(kU)Pi 0) 



W-(kU)R 



(5) 



where P% denotes the probability that Ml is empty 
(m = 0) or occupied (m = 1), while k of its neigh- 
bors are occupied. N is the number of nearest neighbor 
moleucles, [i r is the chemical potential of electrode r, 
U is the Coulomb repulsion between electrons occupying 
nearest neighbor molecules, and the rates are denned by 
W±{E) = r r /±((S + e - fx r )/k B T r ), with f+{x) being 
the fcrmi function and f~(x) = 1 — f + (.r). c is the en- 
ergy of the single spinlcss molecular orbital we take into 
account and r r is the tunnel coupling between Ml and 
electrode r. Since it has no impact on the structure of 
the general theory, we allow for different temperatures in 
the different electrodes, T r ^ T r >. 

Because of the interactions, the current through Ml 
depends on the occupation of its neighbors, i.e., on the 
noncquilibrium charging of the monolayer. Thus, in or- 
der to calculate the current using Eq. , we have to de- 
termine the occupation probability distributions, P^ m ^ ■ 
To do so, we now derive a new type of master (or rate) 
equation. Just as in a normal master equation, the time 
derivative of an occupation probability is given by the 
sum of all tunnel processes enhancing that occupation, 
minus all tunnel processes decreasing that occupation, 
each process being weighted by the occupation probabil- 
ity of the corresponding initial state. We obtain 
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where W A± (£;) = J2,- W ^( E )- The nrst line in E q s - © 
and (|7|) is just like in a standard master equation, de- 
scribing the increase (decrease) in occupation of Ml due 
to tunneling into (out of) it. The second (third and forth) 
lines describe the decrease (increase) in P k due to tun- 
neling into or out of one of M2 (nearest neighbors of Ml), 
which changes k. Here for example P^\, denotes the 
conditional probability that an M2 is occupied and has 
k' occupied neighbors, given that Ml is occupied with k 

I 



occupied nearest neighbors. Note that k' runs from 1 to 
N in the equation for P^, but from to N — 1 in the 
equation for pjf^ . The reason is that if Ml is occupied, 
M2 has to have at least one occupied nearest neighbor, 
and if Ml is unoccupied, M2 can at the most have N — l 
occupied neighbors. 

In the steady-state limit, which is assumed in all results 
presented here and in the main paper, we can set all time 
derivatives to zero in Eqs. ([6]) and (0, as well as in all 
master equations appearing below. 

Equations © and ([7]) have to be supplemented with a 
condition expressing probability normalization: 



3(0) 



^ 1) ) = 1. 



(8) 



The average charge on Ml is given by (q) = —e ^2 k P k ■ 

Clearly, we also have to find the quantities Puul n \ 
which satisfy similar master equations 
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where 7V_ 



N 



1. P/iVl"! for example, denotes the conditional probability that a next-nearest neighbor 
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molecule (M3) is occupied (I = 1) or empty (I — 0), with 
k" occupied nearest neighbors, given that Ml is occupied 
with k occupied nearest neighbors and M2 is occupied 
with k' occupied nearest neighbors. Clearly, a master 
equation for Pj^f]^) will involve even higher order con- 
ditional probabilities, expressing the fact that transport 
through the whole monolayer is in principle correlated. 
Below we will describe approximations which make the 
system of equations close. 

Just as the normal occupation probabilities, the condi- 
tional probabilities must be normalized. This gives rise 
to the additional equations: 



A-l 



1 = 



0) , p (0,l) N 
r k.k' I ' 



E(* 

fc'=0 

JV 



(13) 
(14) 



fc'=i 



In addition, we have 







p(o,i) _ p(i-j) _ p (t,o) _ p (i,i) 

k N — r h. n — * AT l./ — -*n 



fc.O 



N,k' 



O.k' 



(15) 



since, for example, M2 must have at least one occupied 
neighbor if Ml is occupied. 

There is a further possible simplification due to the 
homogeneous nature of the monolayer: Since all M2 
molecules are equivalent, if Ml has k occupied neigh- 
bors, each M2 should be occupied with probability k/N. 
This gives the conditions 



A' 



N 

k 
N 



\ ^ p(m,l) 
/ j r k.k' ' 



fc'=0 
N 



= E n m « 



(m,0) 



(16) 



(17) 



k'=0 



When P^; n) fulfill Eqs. (pj) and ([T7|l. they also auto- 
matically fulfill the standard normalization conditions, 
Eqs. (JTHJl and ([li]) . However, Eqs. (HH) and (JTTJ) give 
twice as many conditions as Eqs. (|13|) and (fl4| . and 
should be used instead of the first line in Eqs. ([§])- (|T2"j). 
corresponding to the terms changing n in Pj^.f 1 " 1 - Thus, 
the occupation of M2 should not be determined by the 
master equation, but is instead fixed by k. We there- 
fore use the master equation only to solve for the k'- 
dependence of P^f 1 ^ ■ 



Mixed mean-field approximation 

We start with a very simple approximation to close 
the system of master equations involving ever higher or- 
ders of conditional probabilities. This approximation is 
not used in any of the results in the main paper, but 
is included here since it provides a simple and intuitive 
method, which is valid for moderate values of the inter- 
action strength, U < ksT. We treat the interactions be- 
tween Ml and M2 exactly, while the interactions between 
M2 and M3 are treated within a mean-field approxima- 
tion. In the master equations for Pj. m \ we should then 
replace k'U by N^(n)U in Eq. (for m = 0) and by 
(1 + N-{n))U in Eq. ((6]) (for m = 1), where (n) is the 
average occupation of the N + 1 molecules Ml and M2, 
given by 



iV 



Nk' ( P, 



(0) 



(18) 



After this approximation, W ± in Eqs. ((6]) and no 
longer depend on k' and we can use Eqs. (fT6| and (fl7|) 
to carry out the sums over the conditional probabilities 

p k"k'' n) - Tllus ' we do not need t0 solvc Ec l s - ©~CH>> 
but can directly solve Eqs. ((6]) and (0, which have now 

become 



P 



(i) 



W + (kU)pjf } - W~{kU)p£ ) 

n[^W- [(1 + N_(n))U] + (l - ^) W+ [(1 + N-(n))U]} if > 

N {^JT W ~ [(1 + X-WW H+i + (l " ^) W + [(1 + iV-(n))C7] P&} , 



(19) 



(0) 



N 



(N-(n)U) + (l-±) 



P 



(0) 



JV 



/c - 1 

N 



W + (N.in^P^ 



(20) 



Together with probability normalization, Eq. ©, this be solved self-consistentlv together with Eq. (|18|) because 
makes a closed set of equations, which, however, has to 
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of the dependence on (n). The resulting occupations, 

Pj, m \ can then be inserted into Eq. ([5]) to calculate the 
current. 



Note that in Eqs. (O and ([20]). and in Eq. © for 
the current, kll appears explicitly in the terms involving 
tunneling into and out of Ml and, in addition, the rele- 
vant interaction for tunneling into M2 is (1 + N-(n))U 
if Ml is occupied [Eq. ([15)1]. but N^(n)U if Ml is unoc- 
cupied [Eq. (JT5J)]. This reflects the fact that the M1-M2 
interactions are treated exactly. Therefore, these master 
equations are capable of producing steps in the noncqui- 
librium monolayer charging and current, for example as 
a function of increasing bias voltage V, and is applicable 
for U < ksT, see Fig. [5] This is in contrast to the much 
simpler standard mean-field approximation, which would 
consist of simply solving only the first line in Eqs. (fT5[) 

and (H3), with kU -> PWU = J2k p k )u i in which case 
the whole equation can be summed over k. The standard 
mean- field approach fails completely unless U <C ksT 
(not shown in Fig. [5]). 

Also note that since we treat the Ml and M2 molecules 
within a different level of approximation, (n) ^ P^\ al- 
though we numerically always find them to be compara- 



ble, unless U ksT. Similarly, if we instead of using 
Eq. ([5]) to calculate the current, calculated the average 
current through Ml and M2, we would obtain a slightly 
different result. 



Truncating the correlations 

We now introduce the approximation scheme which is 
used for the calculations presented in Figs. 2 and 3 in the 
main paper, which does not rely on a mean-field treat- 
ment at any level and is applicable also in the regime 
U > ksT. This is obtained by only taking conditional 
probabilities up to some given order explicitly into ac- 
count. Here, we simply neglect explicit three-charge cor- 
relation terms when solving the set of master equations, 



which means that we replace P^.f'yi — > P^k" (note 



that since we are dealing with conditional probabilities, 
we should not factorize the expressions, i.e., we should 
not include a factor Pu on the right hand side). Now 
the equations for the conditional probabilities, Eqs. (J5J- 
(IT2"1) . become 



J 



k,k' 



(i,i) 



E [w-(k"U)P^l + W + {k"U)P^l 

e"=l 

E [w-{k"u)Pi]^tlw + W+(k"U)P&l lP ^l k ,, 



fe"=l 



(21) 



*& 0) = -N- 



N-l 



£ [w-(k"U)P^l + W+(k"U)Pj ! ?$ 



p(0,0) 
r k,k' 



fc"=0 
N-l 



+ A_ £ [^-(fc^)^?i 1 ^ fe ,+^ + (fc'' f /)^ri^, fc " 



(22) 



N-l 



P&? = -N- 



P (i,o) 

r kM 



53 [w-(k"u)p^l + w+(k"u)p^l 

fc"=0 
N-l 

E [w-(k''u)pi^ +1 p^ + W+(k"U)P^l lP ^ k „ 



(23) 



fe"=0 



6(0,1) 
r k,k< 



-N- 



P 



A_ 



E [w-(k"u)p£il + w + (k"u)pl}$ 

s"=l 

f 

y: [w-{k»u)p^\ x pi^ +w+(k»u)p^i 1 p^ kll 



fe"=l 



(24) 



Together with Eqs. (|16[) and (|17j) . these equations can then inserted into the master equations for the occupa- 
bc solved for Pj^} ■ The conditional probabilities are 
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tion probabilities, Eqs. © and (0, which can be solved 

for Pfr^ , from which we then calculate the current using 
Eq. ©. 

One complication is that Eqs. (f2~T j) - (j!Mj) are non-linear, 

involving terms like Puk' Py^t''- Numerically we can 
deal with this simply by picking some initial distribu- 
tion for the conditional probabilities which originate from 
the truncated higher order terms (Pf?'l),), and then solve 
the equations iteratively, updating this distribution each 
time. Note that what remains are only equations for each 
distribution as a function of fc', i.e., we only need to solve 
for the fc'-dependence of P^k' ■ Thus, we only ever have 
to solve equations of size NxN, although there are 4xiV 
such equations (one for each m, n, and k) which all have 
to be solved within each iteration. As demonstrated in 
Fig. [5j the results are reliable for much larger values of 
U/ksT compared with the mixed mean- field approxima- 
tion introduced above. 

Exact diagonalization of a finite monolayer 

The master equation method described above is suit- 
able only for an infinite homogeneous monolayer. To 
study an inhomogeneous monolayer, e.g., involving de- 
fect sites, we rely instead on a standard master equation 
solution for a finite monolayer. We take a finite num- 
ber of molecules (3x3 molecules on a quadratic lattice 
is used in all results) and diagonalizc exactly the corre- 
sponding monolayer Hamiltonian, Ei -^ml m ^1- (-0 m 
the main paper. Since the total charge in the monolayer, 
Nc, commutes with X^^mlj we can label the many- 
body eigenstates with Nq, and, since we do not consider 
inter-molecular tunneling, each eigenstate can be chosen 
to correspond to a given distribution of these charges 
over the different molecules. We thus label the eigen- 
states \aNc), where a labels the charge configuration. 
The master equation is now solved for the occupation 
probabilities, p a N c i °f these eigenstates, which can then 
be used to calculate the steady-state current from elec- 
trode r into the monolayer 

PaN c = X] ( K K aNc,a'(Nc+s)Pa'(N c +s) 
s=±l a' 

~ K a>(N c +s)MNc PaN c) ' ( 25 ) 

l=^PaJV c , (26) 

aN c 

J r = -e]T Y^ SK a' S (N c +s),aN c PaN C n (27) 

aNc s=±l a' 

where the rate matrix is given by 

^aNca'iNc+s) ~ 1 aN c ,a'{N c +s) 

X I" [(Ea'(N c +s) - E a N c ~ Hr)/kBT r ] , 

(28) 



and KS aN c ,a'(N c +s) = Er K aNc,a'(N c + S y SlnCe We d ° 

not assume a homogeneous monolayer, the tunnel rates, 
^aN c ,a'(N c +s)i depend on both the initial and final many- 
body eigenstates. In particular, T a ^ c a i^ c+S j = when- 
ever \aNc) and \a'(Nc + s)) do not simply differ by the 
occupation of a single molecule, since this is all that can 
be changed by a single SET process. 

Equation (|25j) simply expresses the fact that the 
change in occupation of state \aNc) is determined by 
the sum of all tunnel processes filling \aNc), weighted by 
the occupations of the corresponding initial states, minus 
all tunnel processes emptying \aNc}- In the steady state 
limit, the left hand side of Eq. (|25|) is set to zero. Proba- 
bility normalization is enforced by Eq. (|26[) . The current 
from electrode r, Eq. ([27]) . is given by the sum of all pro- 
cesses involving an electron tunneling into the monolayer 
from electrode r, minus the sum of all processes where 
an electron tunnels out of the monolayer into electrode r, 
weighted by the corresponding occupation probabilities. 

The method described here is used in the final part 
of the main paper to investigate transport in an inho- 
mogeneous molecular monolayer, where some molecules 
couple only weakly to one of the electrodes (and possi- 
bly also have a different LUMO energy). In Fig. [5j we 
use it instead to calculate transport through a homoge- 
neous monolayer, allowing comparison with the infinite 
monolayer methods described above. 



Comparison of different approximations 



In Fig. [5j we compare the result of the different meth- 
ods discussed above for different strengths of the inter- 
molecular Coulomb interaction, U . For U = (not 
shown), all methods give completely identical results. 
Over the whole range of interaction strengths, we find ex- 
cellent agreement between the results obtained by trun- 
cating the correlation, and those from the exact diagonal- 
ization of a finite monolayer. These are the approxima- 
tion schemes used in the main paper. For U ^ fcgT, the 
mixed mean-field solution is clearly inaccurate, predict- 
ing the correct height of the high- voltage current plateau, 
as well as the correct position of the current steps, but 
failing to capture their height, and even erroneously pre- 
dicting NDR for U = 6k B T. Note that this method is 
not used in any of the results in the main paper. 
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FIG. 5. Current (upper panel) and differential conductance (lower panel) as a function of applied bias voltage. In each 
subfigure, results using the three different methods described above are shown: Mixed mean-field approximation (MMF), 
truncation of correlations (CORR) and exact diagonalization of a finite monolayer (ED). We have everywhere used e — 10k B T 
and N = 4 (square lattice). In ED, we used a square lattice with 3x3 molecules, all with the same properties, and periodic 
boundary conditions [using instead open boundary conditions (not shown) does not qualitatively alter the results, but merely 
leads to a small increase in the height of the first peak and reduced height of the satellite peaks, reflecting the fact that 
interactions are somewhat less important in this case]. The strength of the inter-molecular Coulomb interaction is increased, 
with U = k B T in (a), U = 2k B T in (b), U = 4fc s T in (c), and U = Qk B T in (d). 



